rm(list = ls())
library(dplyr)
library(sp)
library(sf)
library(raster)
library(tmap)
library(ggplot2)
library(showtext)
font_add("KaiTi","KaiTi.ttf")
showtext_auto(enable = TRUE)
####################################################################################
setwd("Replication")
####################################################################################

####################################################################################
# Figure 2
####################################################################################
# Figure 2a

map_0 <- read_sf("data/gadm36_MMR_shp/gadm36_MMR_0.shp") %>% 
        st_as_sf(crs = 4326)
map_2 <- read_sf("data/gadm36_MMR_shp/gadm36_MMR_2.shp") %>% 
  st_as_sf(crs = 4326)

map_grid <- st_as_sf(map_0) %>% st_transform(4326)
# create 55km grid - here you can substitute 200 for 50000
grid_55 <- st_make_grid(map_grid, n = c(20, 20),
                        crs = 4326, what = 'polygons') %>% 
            st_intersection(map_grid) %>%
            st_sf(grid_id = 1:length(.))
# create labels for each grid_id
grid_lab <- st_centroid(grid_55) %>%
            cbind(st_coordinates(.))
tm_shape(map_2) + 
        tm_polygons() + tm_layout(frame = FALSE)
tmap_save(filename = "figs/Figure_2a.png", width = 6, height = 8)
##########################################################################

# Figure 2b

grid_400 <- st_make_grid(map_grid, n = c(20, 20),
                         crs = 4326, what = 'polygons') %>% 
            st_sf() %>% 
            mutate(id = row_number())
grid_inter <- st_intersects(map_grid,grid_400,sparse = FALSE)%>%
        as.data.frame() 
grid_select <- grid_inter[,(grid_inter[1,] == TRUE)] 
num <- names(grid_select) %>% 
      substr(start = 2, stop = 100000L) %>% as.numeric()
grid_55 <- grid_400 %>% 
    filter(id %in% num) %>% 
    mutate(id = row_number())

tm_shape(grid_55) + tm_polygons()+
  tm_text("id", size = 0.5)+
  tm_shape(map_0) + tm_borders()

tmap_save(filename = "figs/Figure_2b.png", width = 6, height = 8)

##############################################################################
# Figure 2c
ggplot() +
  geom_sf(data = map_0, fill = 'white', lwd = 0.05) +
  geom_sf(data = grid_55, fill = 'transparent', lwd = 0.3) +
  geom_text(data = grid_lab, aes(x = X, y = Y, label = grid_id), size = 2) +
  coord_sf()+ 
  labs(title = "",
       subtitle = "",
       y = "", x = "")+
  theme(line = element_blank(),rect = element_blank(), #defien the margin line
        axis.text = element_blank(), axis.title = element_blank(), 
        panel.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), axis.ticks.length = unit(0,"cm"),
        panel.spacing = unit(0, "lines"), 
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=16, family = 'KaiTi'),
        plot.margin = unit(c(0,  0, 0, 0), "lines")) 
ggsave("figs/Figure_2b.png", width = 6, height = 8)


#### you can't use this grid for aggregatio because it is only one point rather than the whole area
# read the prio grid data for myanmar
load("data/myanmar_priogrid_cell.RData")
load("data/myanmar_grid.RData")
myanmar_priogrid_cell$grid_id = length(myanmar_priogrid_cell$gid):1
tm_shape(myanmar_priogrid_cell) + tm_polygons()+
    tm_text("grid_id", size = 0.5)+
    tm_shape(map_0) + tm_borders()+
    tm_layout(frame = FALSE)
tmap_save(filename = "figs/Figure_2c.png", width = 6, height = 8)
####################################################################################




####################################################################################
# Figure 3
####################################################################################
# Figure 3a
# light data
light15 <- raster("data/Harmonized_DN_NTL_2015_simVIIRS.tif")
pop15 <- raster("data/gpw_v4_population_count_rev11_2015_30_sec.tif")
# read the prio grid data for myanmar
load("data/myanmar_priogrid_cell.RData")
## convert to SpatialPolygonsDataFrame for overlay
grid <- as_Spatial(myanmar_priogrid_cell)
#make sure they have the same CRS 
proj4string(light15) <- proj4string(grid)
proj4string(pop15) <- proj4string(grid)
# crop the light data by the GRID extent: this will be a squared data
myanmer_light15_map <- crop(light15, extent(grid))
myanmer_pop15_map <- crop(pop15, extent(grid))
# now use the mask function to further crop by the border of the grid
rr <- mask(myanmer_light15_map, grid)
rr2 <- mask(myanmer_pop15_map, grid)
## sum up the total nights by grids
d <- raster::extract(rr, grid, fun = sum, na.rm = TRUE, sp = T )
d2 <- raster::extract(rr2, grid, fun = sum, na.rm = TRUE, sp = T )
# get the aggregated nighttime light data by grid (n = 265)
light15_agg <-  d@data
pop15_agg <-  d2@data
############################################
p_light <- light15_agg %>% 
    dplyr::select(gid, Harmonized_DN_NTL_2015_simVIIRS)
p_pop <- pop15_agg %>% 
  dplyr::select(gid, gpw_v4_population_count_rev11_2015_30_sec)
df_map <- left_join(myanmar_priogrid_cell, p_light, by = "gid")
df_map2 <- left_join(myanmar_priogrid_cell, p_pop, by = "gid")
df_map3 <- left_join(df_map, p_pop, by = "gid")
df_map3 <- df_map3 %>% 
          mutate(gdp_pc = Harmonized_DN_NTL_2015_simVIIRS/gpw_v4_population_count_rev11_2015_30_sec)


p <-  tm_shape(map_0) + tm_borders()+
        tm_shape(df_map) +  
        tm_polygons(col = "Harmonized_DN_NTL_2015_simVIIRS", style = "quantile",
                    title ="2015年夜间灯光") +
        tm_layout(frame = FALSE)

p2 <- tm_shape(map_0) + tm_borders()+
        tm_shape(df_map2) +  
        tm_polygons(col = "gpw_v4_population_count_rev11_2015_30_sec", style = "quantile",
                    title ="2015年人口总数") +
        tm_layout(frame = FALSE)

p3 <- tm_shape(map_0) + tm_borders()+
        tm_shape(df_map3) +  
        tm_polygons(col = "gdp_pc", style = "quantile",
                    title ="2015年人均夜间灯光") +
        tm_layout(frame = FALSE)
tmap_save(tm = p, filename = "figs/Figure_3a.png",width = 6, height = 8)
tmap_save(tm = p2, filename = "figs/Figure_3b.png",width = 6, height = 8)
tmap_save(tm = p3, filename = "figs/Figure_3c.png",width = 6, height = 8)


####################################################################################
# Figure 4
####################################################################################
# Figure 4b
## conflict
load("data/acled_myanmar.RData")
myanmar_battles <- acled_myanmar %>% 
  filter(event_type == "Battles" | event_type == "Explosions/Remote violence")
myanmar_battles <- st_as_sf(myanmar_battles, coords = c("longitude", "latitude"),
                    remove = FALSE, crs = 4326)

p4 <-tm_shape(map_0) + tm_borders()+
  tm_shape (myanmar_battles)+ tm_symbols(col = "blue",
                                 border.lwd = NA, size = 0.1,
                                 shape = 25)+
  tm_grid(labels.size = 2.5,labels.show = T,ticks = F)+
  tm_ylab("纬度（北纬）",size = 5)+
  tm_xlab("经度（东经）",size = 5)+
  tm_layout(fontfamily = "KaiTi")
tmap_save(tm = p4, filename = "figs/Figure_4b.png",width = 6, height = 8)

## Figure 4a
d <- myanmar_battles %>% 
      dplyr::group_by(event_date) %>% 
      dplyr::summarise(events = n()) %>% 
       dplyr::mutate(date = as.Date(event_date, format = "%d %B %Y")) %>% 
       dplyr::mutate(month = as.integer(format(date, "%m")),
                  day = as.integer(format(date, "%d")),
                  yearmonth =  format(date, "%Y-%B")) 

ggplot(d, aes(x = date, y = events)) + geom_line(colour = gray(1/2)) +
  geom_smooth()+
  scale_x_date(date_breaks = "12 months",date_labels = "%m/%y")  +
  labs(y = "冲突数量",x= "时间") + theme_bw() +
  theme(legend.position="none",text = element_text(size=12,family ='KaiTi')) 
ggsave("figs/Figure_4a.png", width = 6, height = 4)


####################################################################################
# Table 1
####################################################################################

vars <- c('battles_sw' , 'battles_s1w' , 'battles_s2w' , 'battles_dum_lag1', 
          'battles_dum_lag2' , 'battles_dum_lag3' ,  'gdppc' , 'log_pop_density',
          'battles_sw' , 'battles_s1w' , 'battles_s2w' , 'battles_dum_lag1', 'battles_dum_lag2',
          'battles_dum_lag3',
          'forest_gc', 'urban_gc', 'rainseas' , 'cmr_mean' , 'water_gc' , 'agri_gc' , 'bdist1' ,  'capdist',
          'battles_dum_lag1' , 'battles_dum_lag2' , 'battles_dum_lag3',
          'battles_sw' , 'battles_s1w' , 'battles_s2w')
vars <- unique(vars)
vars <- sort(vars)
load("data/paper_imputed.RData")
paper_summary <- paper_imputed[which(names(paper_imputed) %in% c("battles_dum",vars))]
paper_summary <- as.data.frame(paper_summary)
library(stargazer)
stargazer(paper_summary, out = "figs/table1.html", type="html",digits = 4,
          summary.stat = c("n", "max", "min", "sd","mean", "median"))
####################################################################################

